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ABSTRACT 


Two  acceptance/rejection  algorithms  for  generating 
random  variates  from  the  beta  distribution  are  developed. 

The  algorithms  use  piece-wise  linear  and  exponential  majorizing 
functions  coupled  with  a piece-wise  linear  minorizing  function. 

The  algorithms  are  exact  to  within  the  accuracy  of  the  computer 
and  are  valid  for  all  parameter  values  greater  than  one.  Execution 
times  are  relatively  insensitive  to  changes  in  parameter  values 
and  are  faster  than  any  previously  published  by  roughly  50%. 
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Generation  of  beta  variates  having  density  function 


f^(x)  - xP  1 (l-x)q_1  / 6(p,q)  0 < x < 1,  l<p,  l<q 

where 

3(p»q)  c /q  up  1 (l-u)q  1 du 

is  considered  here.  Schmeiser  and  Shalaby  [6]  give  three  beta 
algorithms;  B2P,  B4P,  and  BNM;  the  latter  being  a modification 
of  algorithm  BN  of  Ahrens  and  Dieter  [1].  The  computational 
results  in  [6]  indicate  that  in  terms  of  execution  time  these 
three  algorithms  and  algorithm  BB  of  Cheng  [3]  are  each  fastest 
for  some  set  of  parameter  values  (p,q).  In  particular,  BNM  is 
fastest  when  both  p and  q are  large  and  approximately  equal 
(the  distribution  is  approximately  normal),  BB  is  fastest  when 
either  (exclusively)  p or  q is  large  (the  distribution  is  close 
to  gamma),  B2P  is  fastest  when  both  p and  q are  less  than  two, 
and  B4P  is  fastest  otherwise. 

In  this  paper  algorithms  B2P  and  B4P  are  modified  to  yield 
algorithms  B2PE  and  B4PE.  These  algorithms  dominate  B2P  and  B4P, 
respectively,  and  both  dominate  BNM  and  BB  for  all  parameter 
values.  The  algorithms  B2PE  and  B4PE  are  developed  in  Section  1. 
Section  2 contains  computational  results  concerning  execution 
time  and  memory  requirements. 
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1.  ALGORITHM  DEVELOPMENT 


Algorithms  B2P  and  BAP,  as  developed  in  [6],  use 

piece-wise  linear  majorizing  and  minorizing  functions  in  an 

acceptance/rejection  algorithm.  The  efficiency  of  these  algorithms 

is  good  when  both  p and  q are  relatively  small,  but  as  either  or 

both  of  p and  q become  large  the  piece-wise  linear  majorizing 

function  fits  fc(x)  less  well  and  the  algorithms  become  inefficient. 

P 

The  modified  algorithms  developed  here,  B2PE  and  BAPE,  differ  in 

the  logic  used  to  generate  variates  from  the  distribution  tails. 

Following  the  logic  given  in  Schmeiser  [7],  the  linear  majorizing 

function  in  the  tail  is  replaced  by  an  exponential  majorizing  function, 

as  illustrated  in  Figures  1 and  2. 

The  parameters  of  the  exponential  majorizing  function  t (x)  are 

» 

chosen  so  that  f^(Xp)  = t(Xg)  anc*  fg(xQ)  = where  xo  is  tlle 

point  at  which  the  tail  of  the  distribution  begins.  For  B2PE, 
x^,  is  the  point  of  inflection.  For  BAPE,  x^  is  the  point  at  which 

the  line  tangent  to  f-,(x)  at  the  point  of  inflection  intercepts  the 

P 

X axis,  as  shown  in  Figures  1 and  2.  The  validity  of  such  a 
majorizing  function  follows  from  the  Theorem,  which  is  stated  for 
the  right  tail. 

Theorem.  For  any  x^  £ (0,1),  any  K > 0,  and  X « q/d-x^)  - p/x^. 


K Xq  (l-xQ)q  exp [ X ( Xq-x)  ] _>  K xP  (l-x)q  for  all  x E [Xq,  1]. 
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Proof : The  inequality  is  equivalent  to 


xP  (l-x0)q  / exp (-Axq)  _>  xP  (l-x)q  / exp (-Ax)  for  all  x c [xQ,  1]. 


Since  at  x * x^  the  terms  are  equal,  we  need  only  to  show  that 


xP  (l-x)q  exp(Xx)  is  a strictly  decreasing  function  of  x.  The  first 
derivative  with  respect  to  x is 


xP  (l-x)q  exp(Xx)  [p/x  + X - q/ (1-x) ] 


which  is  negative  if  and  only  if 


p/x  - q/ (1-x)  < p/xQ  - q/ (1-Xq) . 


This  condition  is  true  if  and  only  if  p/x  - q/(l-x)  is  a strictly 
decreasing  function  of  x.  Taking  the  first  derivative  yields 


I(l-x)~p  + x2q]  / x2(l-x)2 


which  is  negative,  proving  the  theorem. 


A similar  theorem  holds  for  the  left  tail  using  X - p/x^  “ <3-/(1-Xq). 


Having  shown  the  majorizing  functions  to  be  valid,  we  can  now 
state  the  algorithms.  Each  algorithm  is  composed  of  two  parts: 
set-up  and  generation.  In  the  set-up  the  necessary  constants  are 
calculated  from  the  parameter  values  p and  q.  The  set-up  is  performed 
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only  once  for  each  set  of  parameter  values.  The  generation  logic 
is  executed  each  time  a variate  is  to  be  generated. 


ALGORITHM  B2PE 
Set-up 

1.  Set  P * p-1,  Q = q-1,  R = P+Q,  S - R In  R,  x2  - f2  - f^  - 0, 

x^  = P/R,  x^  = 1.  If  R £ 1,  go  to  4. 

2.  Set  D = (PQ/ (R-l) ) "'■^/R.  If  D _>  x^,  go  to  3.  Otherwise 

set  x2  = x^  - D,  A2  = p/x2  - q/(l-x2), 

f2  = exp (P  In  (x2/P)  + Q In  ((l-x0)/Q)  + S) . 

3.  If  x^  + D _>  1,  go  to  4.  Otherwise  set  x^  = x^  + D, 

= Q/(l-x^)  - p/x^,  f^  = exp (P  In  (x^/P)  + Q In  ((l-x^)/Q)  + S) . 

4.  Set  p^  = x^  - x2>  p2  = ^2^2  + ^1’  ^3  = ^4^4  ^2‘ 

Generation 

5.  Sample  u ^ U(0,1).  Set  u = up^.  Sample  v ^ U(0,1). 

6.  If  u > p^,  go  to  7.  Otherwise  set  x = x2  + u. 

If  x < x^  and  v < f2  + (x-x2> (1— f 2) / (x^-x2 ) , deliver  x. 

If  x _>  x^  and  v < f^  + (x^-x) (1-f^) / (x^-x^) , deliver  x. 

Otherwise  go  to  9. 

7.  If  u > p2,  go  to  8.  Otherwise  set  u * (u-p^) / (p2~p1) , 
x = x2  + ln(u)/A2*  If  v < (A2(x-x2)  + l)/u,  deliver  x. 

If  x <_  0,  go  to  5.  Otherwise  set  v = vf2u  and  go  to  9. 

8.  Set  u = (u-p2)/(p3~p2) , x «=  x^  - In  (u)/A^. 

If  v < (A^(x^-x)  + l)/u,  deliver  x.  If  x > 1,  go  to  5. 

Otherwise  set  v = vf.u. 
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9.  Set  A = In  v.  If  A > -(x-x^)  (R  + R) , go  to  5. 

10.  If  A £ P In  (x/P)  + Q In  ((l-x)/Q)  + S,  deliver  x. 
Otherwise  go  to  5. 


Algorithm  B4PE  uses  a better  piece-wise  fit  which  makes 
it  longer  and  faster. 

ALGORITHM  B4PE 
Set-up 

1.  Set  P = p-1,  Q » q-1,  R = P + Q,  S = R In  R, 

x^  * x2  = 0,  x3  = P/R,  x^  * x^  = 1, 

f1  = f2  = fZ}  xf^  = 0.  If  R < 1,  go  to  4. 

1 /2 

2.  Set  D * (PQ/(R-1))  /R.  If  D > go  to  3.  Otherwise 

set  x9  = x3  - D,  x3  = x2  - ((x2(l-x2))/(P-Rx2)) , 

- P/x3  - Q/ (1-x^ , 

f3  = exp (P  ln(Xj/P)  + Q ln((l-x„)/Q)  + S) . 

f2  - exp (P  ln(x2/P)  + Q ln((l-x2)/Q)  + S). 

3.  If  Xj  + D > 1,  go  to  4.  Otherwise  set 

x^  « x3  + D,  x5  * x^  - ( (x^Cl-x^) ) /(P-Rx^; , 

X5  = Q/(1-x5)  - P/x5, 

f4  “ exp (P  ln(x^/P)  + Q ln((l-x^)/Q)  + S) , 

f5  - exp  (P  ln(x5/P)  + 0 ln((l-x5)/Q)  + S) 

A.  p3  • f 2 (x 3 “ x2^’  Po  “ f4^x4  " x3^  + 

P3  - f3(x2  ~ xi^  + P2*  P4  “ f 5 (x5  “ x^)  + P3> 

P5  - (l-f2)(x3  - x 2 ) + p^,  p^  - (l-f^Mx^  - x3)  + p5» 

p7  “ (f2  ' fl)(x2  ‘ xl)/2  + p6 ’ p8  " (f4  “ f5)(x5  “ x4)/2  + P7» 
P9  “ fl/Xl  + P8’  P10  " f5/X5  + p9‘ 
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Generation 


5. 

Sample 

u % 

U(0,1). 

Set  u = up1Q.  if  u > p , 

go 

to  10. 

6. 

If  u > 

Pl> 

go  to  7. 

Otherwise  deliver  x = x? 

+ 

u/f 

7. 

If  u > 

p2  ’ 

go  to  8. 

Otherwise  deliver  x = x3 

+ 

(u-Pi)/f4. 

8. 

If  u > 

P3* 

go  to  9. 

Otherwise  deliver  x = x^ 

+ 

(u-p2)/f1. 

9. 

Deliver  x 

= x4  + (u- 

P3)/f5- 

10. 

Sample 

W 'V 

U(0,1). 

If  u > py  go  to  11.  Otherwise  set 

x = x2 

+ w(x3  - x2) . 

If  (u  - p4)/(p5  - p4)  <_ 

w, 

deliver  x 

Otherwise 

set  v = f0 

+ (u  - p4>/(x3  - x 2 ) and 

go 

to  16 . 

11. 

If  u > 

P6  ’ 

go  to  12. 

Otherwise  set  x = x3  + w(x 

4 ~ x3)- 

If  (p^  ~ u)/(p  - Pj.)  _>  w,  deliver  x.  Otherwise  set 
v = + (u  - p5)  / (x^  - x3)  and  go  to  16. 


12.  If  u > p.,  go  to  14. 

o 

Otherwise  sample  w0  ^ U(0,1).  If  w^  > w,  set  w = w,, . 

If  u > p^,  go  to  13.  Otherwise  set  x = x^  + w(x^  - x^) , 
v = f^  + 2w(u-p^)/(x2  ~ x^).  If  v < deliver  x. 

Otherwise  go  to  16. 

13.  Set  x = xc  - w(x..  -x,),  v = fc  + 2w(u  - p^)/(xc  - x, ) . 

5 54  5 754 

If  v £ f^w>  deliver  x.  Otherwise  go  to  16. 

14.  If  u > p , go  to  15.  Otherwise  set  u = (p  - u)/(p  - p ), 

x = Xj  + lnCu)/!^.  If  x £ 0,  go  tc  5. 

If  w < (X^(x  - x^)  +l)/u,  deliver  x. 

Otherwise  set  v = wf^u  and  go  to  16. 

15.  Set  u = (p1()  - u)/(p1Q  - pg),  x = x5  - ln(u)/X  . 

If  x _>  1,  go  to  5.  If  w < (X5(x5  - x)  + l)/u,  deliver  x. 
Otherwise  set  v = wf^u. 
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16.  Set  A = In  v.  If  A > - (x  - x^)~(R  + R) , go  to  5. 

17.  If  A <_  P ln(x/P)  + Q ln((l-x)/Q)  + S,  deliver  x. 

Otherwise  go  to  5. 

As  can  be  seen  by  examining  the  algorithms,  composition  is  used 
to  generate  the  variates  from  the  density  corresponding  to  the 
majorizing  function.  In  B2PE  three  regions  compose  the  density, 
while  in  B4PE  ten  regions  are  used.  Their  areas  are  calculated  in 
step  4 in  both  algorithms.  Four  regions  in  B4PE  lie  entirely  under 
fp(x)  and  require  no  rejection  logic.  Thus  in  many  cases  no 
exponential  operations  are  required  to  generate  tb'  beta  variates. 

Step  S of  B2PE  and  step  16  of  B4PE  also  deserve  comment.  Here 
a preliminary  rejection  comparison  is  made  against  the  normal 
majorizing  function  used  by  Ahrens  and  Dieter  [1].  Thus  these 
algorithms  use  both  easy  rejection  and  easy  acceptance  functions. 


2.  COMPUTATIONAL  RESULTS 

The  two  algorithms  developed  in  Section  1 are  compared  to 
algorithms  B2P,  B4P,  and  BB  in  this  section.  Other  beta  variate 
algorithms;  such  as  given  in  Atkinson  and  Pearce  [2],  Fox  [4], 

Jonnk  [5]  and  Whittaker  [8];  are  not  considered  here  since 
it  was  shown  in  [6]  that  the  algorithms  considered  here  dominate 
all  others  in  terms  of  execution  speed. 

The  computational  results  of  this  section  are  based  on  FORTRAN 
implementations  using  the  FTN  compiler  on  Southern  Methodist 
University's  CDC  CYBER  72  computer.  The  U(0,1)  psuedo-random 
numbers  were  generated  by  the  function  subprogram  RANF  intrinsic  in 
the  FTN  compiler. 

The  Table  summarizes  the  computational  results.  Both  B2PE  and 
B4PE  dominate  B2P  and  B4P,  respectively.  In  addition,  both  B2PE 
and  B4PE  dominate  BB,  with  B4PE  being  about  twice  as  fast  as  BB. 

The  Table  does  not  include  timing  information  regarding 
set-up.  For  all  parameter  values  BB  requires  the  least  set-up 
time.  Thus  the  fastest  algorithm  is  dependent  upon  the  number  of 
variates,  M,  to  be  generated  for  fixed  parameter  values.  B4PE 
requires  the  least  total  time  if  any  of  the  follou’ing  three 
conditions  are  satisfied:  1)  M >_  10,  2)  M _>  5 and  Min  (p,  q)  £ 2, 
or  3)  p _<  2 and  q <_  2.  Otherwise  BB  is  fastest. 

In  terms  of  memory  requirements,  B2P  and  B2PE  are  nearly 
equivalent,  as  are  B4P  and  B4PE.  BB  requires  both  the  least  memory 
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TABLE  I 

Comparison  of  the  Beta  Generation  Algorithms 
(in  microseconds) 

Schmeiser  and  Schmeiser  and 


Parameter  Values 

Shalaby 

Babu 

Cheng 

p q 

B2P 

B4P 

B2PE 

B4PE 

BB 

1.0001  1.0001 

.37 

.37 

.37 

.36 

.56 

1.01  1.01 

.36 

.37 

.37 

.37 

.57 

1.5 

.35 

.34 

.36 

.43 

.63 

2. 

.48 

. 46 

.43 

.46 

.65 

5. 

.88 

.42 

.57 

.31 

.72 

10. 

1.7 

.72 

.54 

.32 

.75 

100. 

15.7 

5.7 

.51 

.33 

.79 

1.2  1.2 

.38 

.39 

.38 

.38 

.58 

1.5 

.41 

.41 

. 38 

.39 

.57 

2. 

.41 

.41 

.39 

.41 

.58 

5. 

.45 

.30 

.39 

.28 

.62 

10. 

.78 

.41 

.43 

.30 

.64 

100. 

6.8 

2.1 

.47 

.31 

.70 

2.  2. 

.43 

.44 

.42 

.44 

.61 

5. 

.42 

.33 

.39 

.30 

.63 

10. 

.66 

.38 

.42 

.30 

.67 

100. 

4.8 

1.6 

.48 

.33 

.67 

5.  5. 

.33 

.24 

.33 

.23 

.62 

10. 

.43 

.26 

.35 

.24 

.60 

100. 

2.5 

.82 

.41 

.27 

.59 

10.  10. 

.42 

.27 

.35 

.23 

.65 

100. 

1.8 

.62 

.43 

.25 

.59 

100.  100. 

1.2 

.43 

.35 

.24 

.65 

1000.  1000. 

3.5 

.93 

.36 

.23 

.65 

Memory  Requirements 

498 

655 

467 

690 

335 
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and  Che  least  number  of  lines  of  code. 


No  comparison  of  accuracy  is  made  since  all  algorithms  discussed 
in  this  paper  are  exact  to  within  the  accuracy  of  the  computer. 
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